# Load data
d = read.csv((here("final_analyses", "output", "data", # file path
"respiratory-data.csv")), # file name
header = T)
d = d %>%
as_tibble() %>%
select(!starts_with("X")) %>%
# To create columns with the difference between sample size and events per arm
# This will be useful while using Beta distributions,
# diff = the Beta (b) parameter in a Beta distribution (a,b)
mutate(diff_toci = total_toci - events_toci,
diff_control = total_control - events_control)
Let’s check again how many trials (other than RECOVERY) we have for each outcome and subgroup
d %>%
filter(!str_detect(oxygen, "pooled"),
trial != "RECOVERY") %>%
group_by(outcome, oxygen) %>%
arrange(oxygen) %>%
count() %>%
rename("number of trials" = n) %>%
flextable() %>%
autofit()
outcome | oxygen | number of trials |
discharge | invasive ventilation | 1 |
discharge | non-invasive ventilation | 2 |
discharge | simple oxygen | 1 |
mortality | invasive ventilation | 2 |
mortality | non-invasive ventilation | 3 |
mortality | simple oxygen | 1 |
### To create separate objects per subgroup
# Here I rename columns so these variables are usable in my beta_fun()
# Isolating RECOVERY
recovery_simple_oxygen =
d %>%
filter(trial == "RECOVERY",
oxygen == "simple oxygen") %>%
# Only RECOVERY, so there is no need to sum, just to rename the columns
rename(sum_events_toci = "events_toci",
sum_total_toci = "total_toci",
sum_diff_toci = "diff_toci",
sum_events_control = "events_control",
sum_total_control = "total_control",
sum_diff_control = "diff_control")
recovery_non_invasive =
d %>%
filter(trial == "RECOVERY",
oxygen == "non-invasive ventilation") %>%
# Only RECOVERY, so there is no need to sum, just to rename the columns
rename(sum_events_toci = "events_toci",
sum_total_toci = "total_toci",
sum_diff_toci = "diff_toci",
sum_events_control = "events_control",
sum_total_control = "total_control",
sum_diff_control = "diff_control")
recovery_invasive =
d %>%
filter(trial == "RECOVERY",
oxygen == "invasive ventilation") %>%
# Only RECOVERY, so there is no need to sum, just to rename the columns
rename(sum_events_toci = "events_toci",
sum_total_toci = "total_toci",
sum_diff_toci = "diff_toci",
sum_events_control = "events_control",
sum_total_control = "total_control",
sum_diff_control = "diff_control")
## Group all other trials per outcome and sum the events/sample size
trials_simple_oxygen =
d %>%
filter(trial != "RECOVERY",
oxygen == "simple oxygen") %>%
# Only COVACTA, so no need to sum, just to rename the columns
rename(sum_events_toci = "events_toci",
sum_total_toci = "total_toci",
sum_diff_toci = "diff_toci",
sum_events_control = "events_control",
sum_total_control = "total_control",
sum_diff_control = "diff_control")
trials_non_invasive =
d %>%
filter(trial != "RECOVERY",
oxygen == "non-invasive ventilation") %>%
group_by(outcome) %>%
summarise(sum_events_toci = sum(events_toci),
sum_total_toci = sum(total_toci),
sum_diff_toci = sum(diff_toci),
sum_events_control = sum(events_control),
sum_total_control = sum(total_control),
sum_diff_control = sum(diff_control))
trials_invasive =
d %>%
filter(trial != "RECOVERY",
oxygen == "invasive ventilation") %>%
group_by(outcome) %>%
summarise(sum_events_toci = sum(events_toci),
sum_total_toci = sum(total_toci),
sum_diff_toci = sum(diff_toci),
sum_events_control = sum(events_control),
sum_total_control = sum(total_control),
sum_diff_control = sum(diff_control))
Regarding this outcome, negative risk differences mean benefit.
### Generate posterior beta distributions using a noninformative beta(1,1) prior
set.seed = 123 # set seed for reproducibility (rbeta() function)
n = 10e4 # To determine sampling size
# Here I use the uninformative_posterior_beta() in which you can find in the
# following path:
# here("final_analyses", "script", "functions", "beta_distributions.R")
beta_recovery_simple_oxygen = uniformative_posterior_beta(
recovery_simple_oxygen, # data
"mortality" # outcome
)
### Plot!
# I will use my own functions to standardize my plots
risk_comparison_plot(
beta_recovery_simple_oxygen, # Data object
### start using quotes
"gray50", # Color code for control group
"#D3A464", # Color code for other group
"\nRisk (%)", # X axis label
"RECOVERY trial", # Title
"Simple oxygen only subgroup\n", # Subtitle
### stop using quotes
10, # X axis inferior limit
30, # X axis superior limit
2 # X axis spacing for ticks
)
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.
xlabel = "\nRisk difference (%)"
xlim_inf = -10
xlim_sup = 5
p1 = risk_difference_plot(
beta_recovery_simple_oxygen, # Data object
## start using quotes
"#7EBAC2", # fill color code
xlabel, # X axis label
### stop using quotes
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
2 # X axis label
)
p2 = cumulative_risk_difference_plot(
beta_recovery_simple_oxygen, # Data object
xlabel, # X axis label (in quotes)
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
1, # Spacing between ticks in X axis
"mortality" # Outcome (within quotes)
)
p1 + p2 + plot_annotation(
title = "RECOVERY trial",
subtitle = "Risk difference between groups on simple oxygen only"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.
### Generate posterior beta distributions using a noninformative beta(1,1) prior
set.seed = 123 # set seed for reproducibility (rbeta() function)
n = 10e4 # To determine sampling size
# Here I use the uninformative_posterior_beta() in which you can find in the
# following path:
# here("final_analyses", "script", "functions", "beta_distributions.R")
beta_recovery_non_invasive = uniformative_posterior_beta(
recovery_non_invasive, # data
"mortality" # outcome
)
### Plot!
# I will be using my own functions to standardize my plots
risk_comparison_plot(
beta_recovery_non_invasive, # Data object
### start using quotes
"gray50", # Color code for control group
"#D3A464", # Color code for other group
"\nRisk (%)", # X axis label
"RECOVERY trial", # Title
"Non-invasive ventilation subgroup\n", # Subtitle
### stop using quotes
28, # X axis inferior limit
48, # X axis superior limit
2 # Spacing between ticks in X axis
)
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.
xlabel = "\nRisk difference (%)"
xlim_inf = -12
xlim_sup = 5
p1 = risk_difference_plot(
beta_recovery_non_invasive, # Data object
## start using quotes
"#7EBAC2", # fill color code
xlabel, # X axis label
### stop using quotes
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
2 # Spacing between ticks in X axis
)
p2 = cumulative_risk_difference_plot(
beta_recovery_non_invasive, # Data object
xlabel, # X axis label (in quotes)
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
1, # Spacing between ticks in X axis
"mortality" # Outcome (within quotes)
)
p1 + p2 + plot_annotation(
title = "RECOVERY trial",
subtitle = "Risk difference between groups on non-invasive ventilation"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.
### Generate posterior beta distributions using a noninformative beta(1,1) prior
set.seed = 123 # set seed for reproducibility (rbeta() function)
n = 10e4 # To determine sampling size
# Here I use the uninformative_posterior_beta() in which you can find in the
# following path:
# here("final_analyses", "script", "functions", "beta_distributions.R")
beta_recovery_invasive = uniformative_posterior_beta(
recovery_invasive, # data
"mortality" # outcome
)
### Plot!
# I will be using my own functions to standardize my plots
risk_comparison_plot(
beta_recovery_invasive, # Data object
### start using quotes
"gray50", # Color code for control group
"#D3A464", # Color code for other group
"\nRisk (%)", # X axis label
"RECOVERY trial", # Title
"Invasive mechanical ventilation subgroup\n", # Subtitle
### stop using quotes
36, # X axis inferior limit
58, # X axis superior limit
2 # Spacing between ticks in X axis
) +
# Extra function to better limits the axis in this case
coord_cartesian(c(36, 58))
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.
xlabel = "\nRisk difference (%)"
xlim_inf = -16
xlim_sup = 12
p1 = risk_difference_plot(
beta_recovery_invasive, # Data object
## start using quotes
"#7EBAC2", # fill color code
xlabel, # X axis label
### stop using quotes
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
4 # Spacing between ticks in X axis
)
p2 = cumulative_risk_difference_plot(
beta_recovery_invasive, # Data object
xlabel, # X axis label (in quotes)
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
2, # Spacing between ticks in X axis
"mortality" # Outcome (within quotes)
)
p1 + p2 + plot_annotation(
title = "RECOVERY trial",
subtitle = "Risk difference between groups on invasive mechanical ventilation"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.
priors_simple_oxygen_raw =
d %>%
filter(
trial != "RECOVERY",
oxygen == "simple oxygen"
)
priors_simple_oxygen =
rma(
# Outcome
measure = "RD", # risk difference
# Tocilizumab group
ai = events_toci,
n1i = total_toci,
# Control group
ci = events_control,
n2i = total_control,
data = priors_simple_oxygen_raw %>% filter(outcome == "mortality"),
slab = paste(trial),
method = "REML"
)
forest(priors_simple_oxygen)
set.seed = 123 # set seed for reproducibility (rnorm() function)
n = 10e4 # sampling size
xlabel = "\nRisk difference (%)"
data_prior_posterior_plot(
recovery_simple_oxygen_normal, # Output from normal_approximation()
# start using quotes
"#364D55", # Color for the prior
"#7EBAC2", # Color for RECOVERY
"#A65041", # Color for the posterior
xlabel, # X axis label
"Posterior distribution: RECOVERY and previous evidence", # Title
"Simple oxygen only subgroup", # Subtitle
# stop using quotes
-10, # X axis inferior limit
20, # X axis superior limit
5 # Spacing between ticks in X axis
)
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.
xlim_inf = -8
xlim_sup = 4
p1 = posterior_difference_plot(
recovery_simple_oxygen_normal, # Data object
## start using quotes
"#A65041", # fill color code
xlabel, # X axis label
### stop using quotes
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
2 # Spacing between ticks in X axis
)
p2 = posterior_cumulative_plot(
recovery_simple_oxygen_normal, # Data object
xlabel, # X axis label (in quotes)
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
1, # Spacing between ticks in X axis
"mortality" # Outcome (within quotes)
)
p1 + p2 + plot_annotation(
title = "Posterior distribution: RECOVERY trial and previous evidence",
subtitle = "Simple oxygen only subgroup"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.
priors_noninvasive_raw =
d %>%
filter(
trial != "RECOVERY",
oxygen == "non-invasive ventilation"
)
priors_noninvasive =
rma(
# Outcome
measure = "RD", # risk difference
# Tocilizumab group
ai = events_toci,
n1i = total_toci,
# Control group
ci = events_control,
n2i = total_control,
data = priors_noninvasive_raw %>% filter(outcome == "mortality"),
slab = paste(trial),
method = "REML"
)
forest(priors_noninvasive)
set.seed = 123 # set seed for reproducibility (rnorm() function)
n = 10e4 # sampling size
data_prior_posterior_plot(
recovery_noninvasive_normal, # Output from normal_approximation()
# start using quotes
"#364D55", # Color for the prior
"#7EBAC2", # Color for RECOVERY
"#A65041", # Color for the posterior
xlabel, # X axis label
"Posterior distribution - RECOVERY and previous evidence", # Title
"Non-invasive ventilation subgroup", # Subtitle
# stop using quotes
-15, # X axis inferior limit
5, # X axis superior limit
5 # Spacing between ticks in X axis
)
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.
xlim_inf = -12
xlim_sup = 4
p1 = posterior_difference_plot(
recovery_noninvasive_normal, # Data object
## start using quotes
"#A65041", # fill color code
xlabel, # X axis label
### stop using quotes
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
2 # Spacing between ticks in X axis
)
p2 = posterior_cumulative_plot(
recovery_noninvasive_normal, # Data object
xlabel, # X axis label (in quotes)
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
1, # Spacing between ticks in X axis
"mortality" # Outcome (within quotes)
)
p1 + p2 + plot_annotation(
title = "Posterior distribution: RECOVERY trial and previous evidence",
subtitle = "Non-invasive ventilation subgroup"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.
priors_invasive_raw =
d %>%
filter(
trial != "RECOVERY",
oxygen == "invasive ventilation"
)
priors_invasive =
rma(
# Outcome
measure = "RD", # risk difference
# Tocilizumab group
ai = events_toci,
n1i = total_toci,
# Control group
ci = events_control,
n2i = total_control,
data = priors_invasive_raw %>% filter(outcome == "mortality"),
slab = paste(trial),
method = "REML"
)
forest(priors_invasive)
set.seed = 123 # set seed for reproducibility (rnorm() function)
n = 10e4 # sampling size
data_prior_posterior_plot(
recovery_invasive_normal, # Output from normal_approximation()
# start using quotes
"#364D55", # Color for the prior
"#7EBAC2", # Color for RECOVERY
"#A65041", # Color for the posterior
xlabel, # X axis label
"Posterior distribution - RECOVERY and previous evidence", # Title
"Invasive mechanical ventilation subgroup", # Subtitle
# stop using quotes
-15, # X axis inferior limit
10, # X axis superior limit
5 # Spacing between ticks in X axis
)
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.
xlabel = "\nRisk difference (%)"
xlim_inf = -12
xlim_sup = 8
p1 = posterior_difference_plot(
recovery_invasive_normal, # Data object
## start using quotes
"#A65041", # fill color code
xlabel, # X axis label
### stop using quotes
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
2 # Spacing between ticks in X axis
)
p2 = posterior_cumulative_plot(
recovery_invasive_normal, # Data object
xlabel, # X axis label (in quotes)
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
2, # Spacing between ticks in X axis
"mortality" # Outcome (within quotes)
)
p1 + p2 + plot_annotation(
title = "Posterior distribution: RECOVERY trial and previous evidence",
subtitle = "Invasive mechanical ventilation subgroup"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.
### Create variables for the summary plot
set.seed = 123 # set seed for reproducibility (rnorm())
n = 10e4 # sampling size
## Non-informative priors
delta_beta_recovery_simple_oxygen =
beta_recovery_simple_oxygen %>%
# Calculate difference between groups
summarise("Simple oxygen only" = toci - control)
delta_beta_recovery_non_invasive =
beta_recovery_non_invasive %>%
summarise("Non-invasive ventilation" = toci - control)
delta_beta_recovery_invasive =
beta_recovery_invasive %>%
summarise("Invasive mechanical ventilation" = toci - control)
# Bind all 3 tibbles together
delta_noninformative =
delta_beta_recovery_simple_oxygen %>%
bind_cols(delta_beta_recovery_non_invasive, delta_beta_recovery_invasive)
## Evidence-based priors
delta_recovery_simple_oxygen_normal =
recovery_simple_oxygen_normal %>%
summarise("Simple oxygen only" = rnorm(n,
mean = post.mean,
sd = post.sd)
)
delta_recovery_noninvasive_normal =
recovery_noninvasive_normal %>%
summarise("Non-invasive ventilation" = rnorm(n,
mean = post.mean,
sd = post.sd)
)
delta_recovery_invasive_normal =
recovery_invasive_normal %>%
summarise("Invasive mechanical ventilation" = rnorm(n,
mean = post.mean,
sd = post.sd)
)
# Bind all 3 tibbles together
delta_evidence =
delta_recovery_simple_oxygen_normal %>%
bind_cols(delta_recovery_noninvasive_normal, delta_recovery_invasive_normal)
# Non-informative plot
# X axis
xlim_inf = -10
xlim_sup = 6
ticks_spacing = 2
# Assure subgroup order
# https://stackoverflow.com/questions/12774210/
# how-do-you-specifically-order-ggplot2-x-axis-instead-of-alphabetical-order
level_order = c("Invasive mechanical ventilation", "Non-invasive ventilation", "Simple oxygen only")
p1 = delta_noninformative %>%
# make it tidy for ggplot
pivot_longer(
cols = c("Simple oxygen only":"Invasive mechanical ventilation"),
names_to = "dist", values_to = "samples"
) %>%
ggplot(aes(x = 100 * samples, y = factor(dist, level = level_order))) +
# https://mjskay.github.io/ggdist/articles/slabinterval.html
stat_halfeye(aes(fill = stat(x < 0)),
.width = c(0.5, 0.95)
) +
scale_fill_manual(values = c("gray80", "#7EBAC2")) +
labs(
x = "",
y = "",
title = "Non-informative priors\n"
) +
scale_x_continuous(breaks = seq(from = xlim_inf, to = xlim_sup, ticks_spacing)) +
scale_y_discrete(expand = c(0, 0.3)) +
theme(
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
panel.background = element_blank(),
panel.grid.major.x = element_line(color = "gray80", size = 0.3),
plot.title.position = "plot",
legend.position = "none",
plot.margin = margin(20, 25, 0, 20)
) +
coord_cartesian(xlim = c(xlim_inf, xlim_sup))
# Evidence-based plot
p2 = delta_evidence %>%
# make it tidy for ggplot
pivot_longer(
cols = c("Simple oxygen only":"Invasive mechanical ventilation"),
names_to = "dist", values_to = "samples"
) %>%
ggplot(aes(x = 100 * samples, y = factor(dist, level = level_order))) +
# https://mjskay.github.io/ggdist/articles/slabinterval.html
stat_halfeye(aes(fill = stat(x < 0)),
.width = c(0.5, 0.95)
) +
scale_fill_manual(values = c("gray80", "#A65041")) +
labs(x = "\nRisk difference (%)",
y = "",
title = "Evidence-based priors\n"
) +
scale_x_continuous(breaks = seq(from = xlim_inf, to = xlim_sup, ticks_spacing)) +
scale_y_discrete(expand = c(0, 0.3)) +
theme(
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
panel.background = element_blank(),
panel.grid.major.x = element_line(color = "gray80", size = 0.3),
legend.position = "none",
plot.title.position = 'plot',
plot.margin = margin(0, 25, 10, 20)
) +
coord_cartesian(xlim = c(xlim_inf, xlim_sup))
p1 / p2 + plot_annotation(title = "Posterior distributions on mortality",
theme = theme(plot.title = element_text(size = 20)))
Interval bars depict 50% and 95% credible intervals.
Regarding this outcome, negative risk differences mean benefit.
### Generate posterior beta distributions using a noninformative beta(1,1) prior
set.seed = 123 # set seed for reproducibility (rbeta() function)
n = 10e4 # To determine sampling size
# Here I use the uninformative_posterior_beta() in which you can find in the
# following path:
# here("final_analyses", "script", "functions", "beta_distributions.R")
beta_recovery_simple_oxygen = uniformative_posterior_beta(
recovery_simple_oxygen, # data
"discharge" # outcome
)
### Plot!
# I will use my own functions to standardize my plots
risk_comparison_plot(
beta_recovery_simple_oxygen, # Data object
### start using quotes
"gray50", # Color code for control group
"#5C87C3", # Color code for other group
"\nRisk (%)", # X axis label
"RECOVERY trial", # Title
"Simple oxygen only subgroup\n", # Subtitle
### stop using quotes
58 , # X axis inferior limit
76, # X axis superior limit
2 # X axis spacing for ticks
)
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.
xlabel = "\nRisk difference (%)"
xlim_inf = -2
xlim_sup = 14
p1 = risk_difference_plot(
beta_recovery_simple_oxygen, # Data object
## start using quotes
"#9D95C6", # fill color code
xlabel, # X axis label
### stop using quotes
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
2 # X axis label
)
p2 = cumulative_risk_difference_plot(
beta_recovery_simple_oxygen, # Data object
xlabel, # X axis label (in quotes)
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
1, # Spacing between ticks in X axis
"discharge" # Outcome (within quotes)
)
p1 + p2 + plot_annotation(
title = "RECOVERY trial",
subtitle = "Risk difference between groups on simple oxygen only"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.
### Generate posterior beta distributions using a noninformative beta(1,1) prior
set.seed = 123 # set seed for reproducibility (rbeta() function)
n = 10e4 # To determine sampling size
# Here I use the uninformative_posterior_beta() in which you can find in the
# following path:
# here("final_analyses", "script", "functions", "beta_distributions.R")
beta_recovery_non_invasive = uniformative_posterior_beta(
recovery_non_invasive, # data
"discharge" # outcome
)
### Plot!
# I will be using my own functions to standardize my plots
risk_comparison_plot(
beta_recovery_non_invasive, # Data object
### start using quotes
"gray50", # Color code for control group
"#5C87C3", # Color code for other group
"\nRisk (%)", # X axis label
"RECOVERY trial", # Title
"Non-invasive ventilation subgroup\n", # Subtitle
### stop using quotes
34, # X axis inferior limit
52, # X axis superior limit
2 # Spacing between ticks in X axis
)
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.
xlabel = "\nRisk difference (%)"
xlim_inf = -2
xlim_sup = 16
p1 = risk_difference_plot(
beta_recovery_non_invasive, # Data object
## start using quotes
"#9D95C6", # fill color code
xlabel, # X axis label
### stop using quotes
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
2 # Spacing between ticks in X axis
)
p2 = cumulative_risk_difference_plot(
beta_recovery_non_invasive, # Data object
xlabel, # X axis label (in quotes)
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
1, # Spacing between ticks in X axis
"discharge" # Outcome (within quotes)
)
p1 + p2 + plot_annotation(
title = "RECOVERY trial",
subtitle = "Risk difference between groups on non-invasive ventilation"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.
### Generate posterior beta distributions using a noninformative beta(1,1) prior
set.seed = 123 # set seed for reproducibility (rbeta() function)
n = 10e4 # To determine sampling size
# Here I use the uninformative_posterior_beta() in which you can find in the
# following path:
# here("final_analyses", "script", "functions", "beta_distributions.R")
beta_recovery_invasive = uniformative_posterior_beta(
recovery_invasive, # data
"discharge" # outcome
)
### Plot!
# I will be using my own functions to standardize my plots
risk_comparison_plot(
beta_recovery_invasive, # Data object
### start using quotes
"gray50", # Color code for control group
"#5C87C3", # Color code for other group
"\nRisk (%)", # X axis label
"RECOVERY trial", # Title
"Invasive mechanical ventilation subgroup\n", # Subtitle
### stop using quotes
10, # X axis inferior limit
26, # X axis superior limit
2 # Spacing between ticks in X axis
) +
# Extra function to better limits the axis in this case
coord_cartesian(c(10, 26))
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.
xlabel = "\nRisk difference (%)"
xlim_inf = -8
xlim_sup = 12
p1 = risk_difference_plot(
beta_recovery_invasive, # Data object
## start using quotes
"#9D95C6", # fill color code
xlabel, # X axis label
### stop using quotes
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
4 # Spacing between ticks in X axis
)
p2 = cumulative_risk_difference_plot(
beta_recovery_invasive, # Data object
xlabel, # X axis label (in quotes)
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
2, # Spacing between ticks in X axis
"discharge" # Outcome (within quotes)
)
p1 + p2 + plot_annotation(
title = "RECOVERY trial",
subtitle = "Risk difference between groups on invasive mechanical ventilation"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.
priors_simple_oxygen_raw =
d %>%
filter(
trial != "RECOVERY",
oxygen == "simple oxygen"
)
priors_simple_oxygen =
rma(
# Outcome
measure = "RD", # risk difference
# Tocilizumab group
ai = events_toci,
n1i = total_toci,
# Control group
ci = events_control,
n2i = total_control,
data = priors_simple_oxygen_raw %>% filter(outcome == "discharge"),
slab = paste(trial),
method = "REML"
)
forest(priors_simple_oxygen)
set.seed = 123 # set seed for reproducibility (rnorm() function)
n = 10e4 # sampling size
xlabel = "\nRisk difference (%)"
data_prior_posterior_plot(
recovery_simple_oxygen_normal, # Output from normal_approximation()
# start using quotes
"#364D55", # Color for the prior
"#9D95C6", # Color for RECOVERY
"#DBA237", # Color for the posterior
xlabel, # X axis label
"Posterior distribution: RECOVERY and previous evidence", # Title
"Simple oxygen only subgroup", # Subtitle
# stop using quotes
-10, # X axis inferior limit
20, # X axis superior limit
5 # Spacing between ticks in X axis
)
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.
xlim_inf = 0
xlim_sup = 14
p1 = posterior_difference_plot(
recovery_simple_oxygen_normal, # Data object
## start using quotes
"#DBA237", # fill color code
xlabel, # X axis label
### stop using quotes
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
2 # Spacing between ticks in X axis
)
p2 = posterior_cumulative_plot(
recovery_simple_oxygen_normal, # Data object
xlabel, # X axis label (in quotes)
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
1, # Spacing between ticks in X axis
"discharge" # Outcome (within quotes)
)
p1 + p2 + plot_annotation(
title = "Posterior distribution: RECOVERY trial and previous evidence",
subtitle = "Simple oxygen only subgroup"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.
priors_noninvasive_raw =
d %>%
filter(
trial != "RECOVERY",
oxygen == "non-invasive ventilation"
)
priors_noninvasive =
rma(
# Outcome
measure = "RD", # risk difference
# Tocilizumab group
ai = events_toci,
n1i = total_toci,
# Control group
ci = events_control,
n2i = total_control,
data = priors_noninvasive_raw %>% filter(outcome == "discharge"),
slab = paste(trial),
method = "REML"
)
forest(priors_noninvasive)
set.seed = 123 # set seed for reproducibility (rnorm() function)
n = 10e4 # sampling size
data_prior_posterior_plot(
recovery_noninvasive_normal, # Output from normal_approximation()
# start using quotes
"#364D55", # Color for the prior
"#9D95C6", # Color for RECOVERY
"#DBA237", # Color for the posterior
xlabel, # X axis label
"Posterior distribution - RECOVERY and previous evidence", # Title
"Non-invasive ventilation subgroup", # Subtitle
# stop using quotes
-10, # X axis inferior limit
15, # X axis superior limit
5 # Spacing between ticks in X axis
)
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.
xlim_inf = -2
xlim_sup = 14
p1 = posterior_difference_plot(
recovery_noninvasive_normal, # Data object
## start using quotes
"#DBA237", # fill color code
xlabel, # X axis label
### stop using quotes
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
2 # Spacing between ticks in X axis
)
p2 = posterior_cumulative_plot(
recovery_noninvasive_normal, # Data object
xlabel, # X axis label (in quotes)
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
1, # Spacing between ticks in X axis
"discharge" # Outcome (within quotes)
)
p1 + p2 + plot_annotation(
title = "Posterior distribution: RECOVERY trial and previous evidence",
subtitle = "Non-invasive ventilation subgroup"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.
priors_invasive_raw =
d %>%
filter(
trial != "RECOVERY",
oxygen == "invasive ventilation"
)
priors_invasive =
rma(
# Outcome
measure = "RD", # risk difference
# Tocilizumab group
ai = events_toci,
n1i = total_toci,
# Control group
ci = events_control,
n2i = total_control,
data = priors_invasive_raw %>% filter(outcome == "discharge"),
slab = paste(trial),
method = "REML"
)
forest(priors_invasive)
set.seed = 123 # set seed for reproducibility (rnorm() function)
n = 10e4 # sampling size
data_prior_posterior_plot(
recovery_invasive_normal, # Output from normal_approximation()
# start using quotes
"#364D55", # Color for the prior
"#9D95C6", # Color for RECOVERY
"#DBA237", # Color for the posterior
xlabel, # X axis label
"Posterior distribution - RECOVERY and previous evidence", # Title
"Invasive mechanical ventilation subgroup", # Subtitle
# stop using quotes
-10, # X axis inferior limit
25, # X axis superior limit
5 # Spacing between ticks in X axis
)
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.
xlabel = "\nRisk difference (%)"
xlim_inf = -8
xlim_sup = 12
p1 = posterior_difference_plot(
recovery_invasive_normal, # Data object
## start using quotes
"#DBA237", # fill color code
xlabel, # X axis label
### stop using quotes
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
4 # Spacing between ticks in X axis
)
p2 = posterior_cumulative_plot(
recovery_invasive_normal, # Data object
xlabel, # X axis label (in quotes)
xlim_inf, # X axis inferior limit
xlim_sup, # X axis superior limit
2, # Spacing between ticks in X axis
"discharge" # Outcome (within quotes)
)
p1 + p2 + plot_annotation(
title = "Posterior distribution: RECOVERY trial and previous evidence",
subtitle = "Invasive mechanical ventilation subgroup"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.
### Create variables for the summary plot
set.seed = 123 # set seed for reproducibility (rnorm())
n = 10e4 # sampling size
## Non-informative priors
delta_beta_recovery_simple_oxygen =
beta_recovery_simple_oxygen %>%
# Calculate difference between groups
summarise("Simple oxygen only" = toci - control)
delta_beta_recovery_non_invasive =
beta_recovery_non_invasive %>%
summarise("Non-invasive ventilation" = toci - control)
delta_beta_recovery_invasive =
beta_recovery_invasive %>%
summarise("Invasive mechanical ventilation" = toci - control)
# Bind all 3 tibbles together
delta_noninformative =
delta_beta_recovery_simple_oxygen %>%
bind_cols(delta_beta_recovery_non_invasive, delta_beta_recovery_invasive)
## Evidence-based priors
delta_recovery_simple_oxygen_normal =
recovery_simple_oxygen_normal %>%
summarise("Simple oxygen only" = rnorm(n,
mean = post.mean,
sd = post.sd)
)
delta_recovery_noninvasive_normal =
recovery_noninvasive_normal %>%
summarise("Non-invasive ventilation" = rnorm(n,
mean = post.mean,
sd = post.sd)
)
delta_recovery_invasive_normal =
recovery_invasive_normal %>%
summarise("Invasive mechanical ventilation" = rnorm(n,
mean = post.mean,
sd = post.sd)
)
# Bind all 3 tibbles together
delta_evidence =
delta_recovery_simple_oxygen_normal %>%
bind_cols(delta_recovery_noninvasive_normal, delta_recovery_invasive_normal)
# Non-informative plot
# X axis
xlim_inf = -8
xlim_sup = 14
ticks_spacing = 2
# Assure subgroup order
# https://stackoverflow.com/questions/12774210/
# how-do-you-specifically-order-ggplot2-x-axis-instead-of-alphabetical-order
level_order = c("Invasive mechanical ventilation", "Non-invasive ventilation", "Simple oxygen only")
p1 = delta_noninformative %>%
# make it tidy for ggplot
pivot_longer(
cols = c("Simple oxygen only":"Invasive mechanical ventilation"),
names_to = "dist", values_to = "samples"
) %>%
ggplot(aes(x = 100 * samples, y = factor(dist, level = level_order))) +
# https://mjskay.github.io/ggdist/articles/slabinterval.html
stat_halfeye(aes(fill = stat(x > 0)),
.width = c(0.5, 0.95)
) +
scale_fill_manual(values = c("gray80", "#9D95C6")) +
labs(
x = "",
y = "",
title = "Non-informative priors\n"
) +
scale_x_continuous(breaks = seq(from = xlim_inf, to = xlim_sup, ticks_spacing)) +
scale_y_discrete(expand = c(0, 0.3)) +
theme(
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
panel.background = element_blank(),
panel.grid.major.x = element_line(color = "gray80", size = 0.3),
plot.title.position = "plot",
legend.position = "none",
plot.margin = margin(20, 25, 0, 20)
) +
coord_cartesian(xlim = c(xlim_inf, xlim_sup))
# Evidence-based plot
p2 = delta_evidence %>%
# make it tidy for ggplot
pivot_longer(
cols = c("Simple oxygen only":"Invasive mechanical ventilation"),
names_to = "dist", values_to = "samples"
) %>%
ggplot(aes(x = 100 * samples, y = factor(dist, level = level_order))) +
# https://mjskay.github.io/ggdist/articles/slabinterval.html
stat_halfeye(aes(fill = stat(x > 0)),
.width = c(0.5, 0.95)
) +
scale_fill_manual(values = c("gray80", "#DBA237")) +
labs(x = "\nRisk difference (%)",
y = "",
title = "Evidence-based priors\n"
) +
scale_x_continuous(breaks = seq(from = xlim_inf, to = xlim_sup, ticks_spacing)) +
scale_y_discrete(expand = c(0, 0.3)) +
theme(
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
panel.background = element_blank(),
panel.grid.major.x = element_line(color = "gray80", size = 0.3),
legend.position = "none",
plot.title.position = 'plot',
plot.margin = margin(0, 25, 10, 20)
) +
coord_cartesian(xlim = c(xlim_inf, xlim_sup))
p1 / p2 + plot_annotation(title = "Posterior distributions on hospital discharge",
theme = theme(plot.title = element_text(size = 20)))
Interval bars depict 50% and 95% credible intervals.
sessionInfo()
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] metafor_2.4-0 Matrix_1.3-2 here_1.0.1 patchwork_1.1.1
## [5] ggdist_2.4.0 tidybayes_2.3.1 flextable_0.6.4 forcats_0.5.1
## [9] stringr_1.4.0 dplyr_1.0.5 purrr_0.3.4 readr_1.4.0
## [13] tidyr_1.1.3 tibble_3.1.1 ggplot2_3.3.3 tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-152 fs_1.5.0 lubridate_1.7.10
## [4] httr_1.4.2 rprojroot_2.0.2 tools_4.0.4
## [7] backports_1.2.1 bslib_0.2.4 utf8_1.2.1
## [10] R6_2.5.0 DBI_1.1.1 colorspace_2.0-0
## [13] withr_2.4.1 tidyselect_1.1.0 curl_4.3
## [16] compiler_4.0.4 cli_2.4.0 rvest_1.0.0
## [19] arrayhelpers_1.1-0 xml2_1.3.2 officer_0.3.17
## [22] sass_0.3.1 scales_1.1.1 systemfonts_1.0.1
## [25] digest_0.6.27 rmarkdown_2.7 base64enc_0.1-3
## [28] pkgconfig_2.0.3 htmltools_0.5.1.1 highr_0.9
## [31] dbplyr_2.1.1 fastmap_1.1.0 rlang_0.4.10
## [34] readxl_1.3.1 rstudioapi_0.13 jquerylib_0.1.3
## [37] generics_0.1.0 farver_2.1.0 svUnit_1.0.3
## [40] jsonlite_1.7.2 zip_2.1.1 distributional_0.2.2
## [43] magrittr_2.0.1 Rcpp_1.0.6 munsell_0.5.0
## [46] fansi_0.4.2 gdtools_0.2.3 lifecycle_1.0.0
## [49] stringi_1.5.3 yaml_2.2.1 plyr_1.8.6
## [52] grid_4.0.4 crayon_1.4.1 lattice_0.20-41
## [55] haven_2.4.0 hms_1.0.0 knitr_1.32
## [58] pillar_1.6.0 uuid_0.1-4 reprex_2.0.0
## [61] glue_1.4.2 evaluate_0.14 V8_3.4.0
## [64] data.table_1.14.0 modelr_0.1.8 vctrs_0.3.7
## [67] cellranger_1.1.0 gtable_0.3.0 assertthat_0.2.1
## [70] cachem_1.0.4 xfun_0.22 broom_0.7.6
## [73] coda_0.19-4 memoise_2.0.0 ellipsis_0.3.1
## [76] rdocsyntax_0.4.1.9000